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Abstract 

This paper proposes a novel computationally efficient dynamic bi-orthogonality based approach for cali- 
bration of a computer simulator with high dimensional parametric and model structure uncertainty. The 
proposed method is based on a decomposition of the solution into mean and a random field using a generic 
Karhunnen-Loeve expansion. The random field is represented as a convolution of separable Hilbert spaces 
in stochastic and spatial dimensions that are spectrally represented using respective orthogonal bases. In 
particular, the present paper investigates generalized Polynomial Chaos bases for stochastic dimension and 
eigenfunction bases for spatial dimension. Dynamic orthogonality is used to derive closed form equations for 
the time evolution of mean, spatial and the stochastic fields. The resultant system of equations consists of a 
partial differential equation (PDE) that define dynamic evolution of the mean, a set of PDEs to define the 
time evolution of eigenfunction bases, while a set of ordinary differential equations (ODEs) define dynamics 
of the stochastic field. This system of dynamic evolution equations efficiently propagates the prior para- 
metric uncertainty to the system response. The resulting bi-orthogonal expansion of the system response 
is used to reformulate the Bayesian inference for efficient exploration of the posterior distribution. Efficacy 
of the proposed method is investigated for calibration of a 2D transient diffusion simulator with uncertain 
source location and diffusivity. Computational efficiency of the method is demonstrated against a Monte 
Carlo method and a generalized Polynomial Chaos approach. 

Keywords: Bayesian Framework, Dynamically Bi-orthogonal Field Equations, Karhunnen-Loeve 
Expansion, Generalized Polynomial Chaos Basis 



1. Introduction 

Recent advancements in digital technologies have facilitated the use of computer simulators for investi- 
gation of large scale systems. However, computer simulators are fraught with uncertainties due to poorly 
known/unknown model, parameters, initial and boundary conditions etc. [15]. Various researchers have in- 
vestigated effect of these uncertainties on the credibility of a computer simulator and established uncertainty 
quantification and calibration as an integral aspect of a modeling and simulation process [28, 29, 32, 27]. 
This paper focuses on the Bayesian approach that provide a formal framework to identify, characterize and 
quantify the uncertainties, and provides a generic inference method for calibration of a computer simulator 
using limited and noisy experimental data [12, 3, 11, 13, 7]. 

The Bayesian framework is preferred over more traditional calibration methods due to its ability to 
provide complete posterior statistics of the parameters of interest. Sampling techniques such as Markov 
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Chain Monte Carlo (MCMC) method [9, 2] are used for exploration of the posterior statistics, especially for 
calibration of nonlinear dynamical simulators in non-Gaussian settings. Satisfactory approximation of the 
posterior distribution and associated statistics using MCMC requires evaluation of the simulator at large 
number of input settings, often in the range of 10 3 — 10 6 . Collection of large number of samples become 
computationally prohibitive for simulation of a large scale system, imposing key challenge for implementation 
of the Bayesian framework. To make the Bayesian framework accessible to large-scale problems, it is 
necessary to develop computationally efficient uncertainty propagation and calibration techniques. 

Marzouk et al. [33] have proposed computationally efficient implementation of the Bayesian framework 
using stochastic spectral methods. Stochastic spectral projection (SSP) based methods are extensively 
used for uncertainty propagation as a computationally efficient alternative to Monte Carlo methods with 
comparable accuracy. Homogeneous Chaos theory introduced by Wiener [16, 17] is an earliest exposition 
of SSP method, where random variables are represented as an expansion series in orthogonal Hermite 
polynomials that converges in mean-square sense [23]. Present state of the art in the field of SSP based 
methods for uncertainty propagation is based on the generalized Polynomial Chaos (gPC) method. The 
method has been successfully implemented for solution of stochastic finite element methods [21, 19, 22] 
and stochastic fluid flow problems [18, 5]. Xiu and Karniadakis [6] have extended the method to a set of 
Askey scheme of orthogonal polynomials. Subsequently, the method has been applied by various researchers 
for uncertainty propagation through simulators of systems of engineering importance [4, 10, 30, 8]. The 
Bayesian calibration formulation proposed by Marzouk et al. [33] uses the gPC method to propagate the 
prior parametric uncertainty to the simulator prediction. Resultant gPC expansion of the prediction is used 
in the Bayes theorem to define the likelihood. The methodology is further extended by Marzouk and Najm 
[34] for inference of spatially /temporally varying uncertain parameters. 

Although the gPC method provides computationally efficient estimation of the uncertainty, computa- 
tional cost of the implementation grows significantly as the number of stochastic dimensions increases [25]. 
Such a high dimensional uncertainty typically arises for a simulator with a large number of uncertain pa- 
rameters, and more predominantly in the case of a spatially /temporally varying uncertain parameter with 
rapidly decaying covariance functions. The research work presented in this paper addresses the Bayesian 
framework for calibration of a simulator with high dimensional uncertainty. 

Sapsis and Lermusiaux [24] have proposed dynamically orthogonal field equations (DOFE) method for 
efficient propagation of high-dimensional uncertainty. The method uses decomposition of the system response 
into a mean and stochastic dynamical component using a truncated generalized Karhunnen-Loeve expansion. 
The stochastic component is spectrally represented in terms of orthogonal eigenfunction basis in spatial 
dimension, while the respective coefficients define the time varying stochastic dimension. The Dynamic 
Orthogonality (DO) condition [24] is used to derive the closed form evolution equations for the mean, 
eigenfunction basis and the stochastic coefficients. 

However, the dynamic evolution equations of Sapsis and Lermusiaux [24] does not impose any geometric 
structure on stochastic dimension, which makes it hard to directly apply the DOFE methodology for Bayesian 
calibration problems. This paper proposes a dynamic &i-orthogonality-based approach that extends the 
DOFE method for the Bayesian calibration of a computer simulator with high dimensional uncertainty. The 
proposed bi-orthogonal method uses spectral expansion of the stochastic field in the gPC basis, imposing 
geometric structure on the stochastic dimension. The random coefficients of the truncated generalized 
Karhunnen-Loeve expansion, obtained using dynamically orthogonal field equations, are projected on the 
gPC basis using Galerkin projection [22]. The resultant field equations are termed here as dynamically bi- 
orthogonal field equations (DBFE). The Bayesian calibration approach proposed in this paper uses DBFE 
to project the prior parametric uncertainty to the system response. The resultant bi-orthogonal expansion 
is used in the likelihood during MCMC sampling for Bayesian inference of the uncertain parameters. The 
proposed DBFE method provides substantial computational speedup over gPC based method for Bayesian 
calibration. Efficacy of the proposed method is demonstrated for calibration of a 2D transient diffusion 
equation simulator with uncertain source location and the diffusivity field. Note that a preliminary version 
of this work was reported in [26] , while this article includes (a) extension of the method to take into account 
model structural uncertainty; (b) substantially refined and expanded theoretical analysis; and (c) additional 
numerical results demonstrating computational efficiency of the proposed methodology. 
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The rest of the paper is organized as follows. In section 2, statistical formulation is discussed in detail. 
The proposed DBFE-based Bayesian method is presented in section 3. Section 4 provides numerical results 
for 2D transient diffusion equation. Finally, the paper is summarized and concluded in section 5. 



2. Statistical Formulation 

The proposed Bayesian framework is developed for a simulator T(x, t, 9(u))), where t € R>o is time, 
x G X C M d is a spatial dimension and 9 (id) is a set of uncertain parameters that induces uncertainty in 
the predictions. Note that the simulator T(x,t,9(u))) is defined over a probability space (f2,.F, "P), where 
uj € O is a set of elementary events, JF is associated er-algebra and V is a probability measure defined over 
J- . In this paper, the proposed method is particularly developed for simulators with a model given by the 
partial differential equation of type 

Mx gt t;Uj) = C [<x, t; W ); 0(w)] , (SPDE) 

where u(x,t;iu) is the system response and C is an arbitrary differential operator. Equation (SPDE) 
is known as stochastic partial differential equation (SPDE). (SPDE) is initialized using a random field 
u(x, 0; w), while, the boundary condition is given by 

B(fi, t; lo) = h(0, t;u); /3 e cU\ w e O, (1) 

where B is a linear differential operator. 

The simulator T(x,t, 9(oj)) approximates the physical phenomena within the limits of available knowl- 
edge. Let C(x, t) be the 'true' but unknown model that perfectly represents the physical phenomena. Since 
T(x, t, 9(oj)) is an approximate representation of the physical phenomena, the simulator predictions deviates 
from ((x,t) by S(x,t), where S(x,t) is known as the discrepancy function. The relationship between ((x,t) 
and T(x, t, 0(u)) is given by [12] 

C(x,t) = T(x,t,9(cj))+S(x,t), (2) 

where 9(u>) denotes the 'true' value of the uncertain parameters. 

The proposed Bayesian framework uses experimental observations at finite locations for inference of 
the uncertain parameters. At time t, let the system be experimentally observed at M spatial locations 
{xf, i = 1, . . . , M}. The measurement at a^, denoted as y e {xi, t), is given by 

y e (xi,t) = ((xi,t) + c(xi,t), (3) 

where e(xi,t) is the measurement uncertainty. Let y e = {y e (xi,t); i = 1, ...,M} be the set of available 
experimental observations. Using y e , the uncertain parameters 9 and the discrepancy function 5(x,t) can 
be inferred through the Bayes theorem as 

/ (§{u),6(x, t) \y e )<xf (y e \ T(x, t, 9(u)),5(x, i)) x / (§(u),5(x, i)) , (4) 

where f(9(uj),5(x,t)) is the prior, f(y e \ T(x,t,9(uj)),S(x,t)) is the likelihood and f(9(uj),S(x,t) \ y e ) is 
the posterior probability distribution. 

Uncertainty in the experimental observations, e(x 1 t), is assumed to be specified using a zero mean 
Gaussian distribution with covariance matrix 

E e = a\l M , (5) 

where / is the M x M identity matrix. In this paper, the proposed Bayesian framework is developed for 
independent prior uncertainty in 9(uj) and 5(x,t) 1 with the prior for 5(x,t) given by a zero-mean Gaussian 
process with a covariance function of the form 

Y lS {x l ,x 2 ) = a} exp I - V \ t (x\ - x l 2 f , (6) 



where a\ is the variance and A^ is the correlation strength of the covariance function, which are treated as 
uncertain hyper-parameters. In the present paper, inverse Gamma prior /G(a CT , j5 a ) is used for er|, while, the 
Gamma prior G(a\ i , (3^) is used for A^ [12, 20, 1]. The uncertainty in S(x,t) is specified using hierarchical 
zero-mean Gaussian process prior as 

f(5(x,t)^l\ t )^\Y, s \-^ exp (-Vl^) x (of)- 3 *- ^xp^-^ x J[(X l ) a ^- 1 exp (-/3 Xi X t ) (7) 

where d = {S(xi,t); i — 1,...,N}. Use (7) and (5) in the Bayes theorem (4) and marginalize 5(x,t) to 
obtain 

f(0(uj),a 2 s ,X l ) oc| S \~i exp (-I^sr 1 *^ x (of)- '" 1 exp (-^) x Y[(X t r^ 1 exp (-/3 Ai A,) x f(§(u)), 
where £ = £5 + £ e and J7 = {y e (:Ci, i) — T(x t , i, 0(uj)); i = 1, A/}. 



3. Proposed Methodology 

Equation (8) can be solved by sampling from the posterior distribution using MCMC, which requires 
evaluation of T{xi,t 1 0{uj)) for each sample, which is computationally prohibitive for large-scale system 
simulators. The approach proposed in this paper requires single evaluation of (SPDE) using dynamically 
bi-orthogonal field equations. The DBFE is used for propagating the prior uncertainty in 0{uj) to the system 
response. The resultant bi-orthogonal expansion of the system response is used in (8) to define the posterior 
distribution, which is explored using the MCMC. The proposed method is described in detail in this section. 

3.1. Dynamically Bi-orthogonal Field Equations 

The proposed DBFE method is based on the dynamically orthogonal held equations (DOEF) proposed 
by Sapsis and Lermusiaux [24]. Consider a generic Karhunnen-Loeve expansion of u{x,t;ui) truncated at 
N terms as 



JV 

n\ ' - 

i=l 

' 2/ 



i(x, t; uj) = u(x, t) Y l (t; u))ui(x, t), (9) 



where u(x,t) is the mean, u.i(x,t) are the functions that form the complete orthonormal basis on L (X), 
while Yi(t; lo) are the zero-mean independent random variables. Note that throughout this paper the equality 
sign, =, is used to represent the approximate equality, if no confusion is expected. Substituting the expansion 
(9) in (SPDE) gives 

du(x,t) / ,dYi(t:uj) v-^-,,, .duAx,t) „. , . „. .. . . 

— ^-L + Ui(x, t) ^ ' > + ]T Yi(t; cj) ^ ; = C[u(x, t; «); 0(u>)]. (10) 

i=l i=l 

Note that the quantities u(x, t), Yi(t; u>) and Ui(x, t) are dependant on each other through (9). Thus, varying 
u(x,t), Yi(t;uj) and Ui(x,t) concurrently makes (10) redundant, necessitating imposition of the additional 
constraint to derive the independent evolution equations for unknown quantities. 

Sapsis and Lermusiaux [24] proposed imposition of dynamic orthogonality (DO) condition to derive the 
independent evolution equations for u(x,t), Yi(t;uj) and Ui(x,t). The DO condition constraints the time 
evolution of u.i(x, t) such that 



dUiM Uj (x,t)} =0: 1i,j = l,...,N, (11) 
X 



dt 



where (•, •) is inner product. Note that the DO condition ensures that Ui(x,t) preserves orthonormality over 
the time evolution of (SPDE). 
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Remark 1. In this paper, inner product is defined over spatial and stochastic dimensions. Inner product 
over the spatial dimension is defined as 

(u(x,t;u)),v(x, t;u)) x = / u(x, t; u)v(x, t; u))dx, (12) 

Jx 

while the inner product over the stochastic dimension is defined as 

{u(x,t;cj),v(x,t;u})) Q = / u(x,t;u>)v(x,t;u)ctP(u)). (13) 

Jo. 

Using the DO condition, the independent evolution equations for u{x, t), Ui{x, t) and Yi(t; u) are derived 
as follows [24]. 

3.1.1. Dynamically Orthogonal Field Equations 

Apply the expectation operator to (10) to obtain the evolution equations for u(x,t) as 

^^=E»[£[u(x,t ] u) ] e(u)}}. (14) 
Multiply (10) by Yj(t;uj) and apply the expectation operator to have 

V cw-m dUt< ?' t] + E c ^w^) Ui ^ t] = E " l£[<^t;");0(")}Yj(t;")} > (is) 

i=l i=l dt 3 

where Cy^Y^i) denote the covariance between Yi(t;ui) and Yj(t;uj). By multiplying Uk(x,t) to (15), taking 
the inner product and applying the expectation operator gives 

C dY M t) Y ] (t) = ^ w I' c ( u ( x '*;' j );' j ) y j(*! a; )] ^ u k{x,t)) x , (i6) 

which on substitution in (15) provides 



N 



dui(x, t) 

J Ydt)Y 3 (t) 



(17) 



/V 



= E u [£[u(a:,t,u;);0(u;)]l$(t;w)] - ^ [C[u(x, t, u); e{u)}Yj(t; u>)} , u k (x, t)) x u k (x, t). 

k=l 

Equation (17) can be written in the matrix form 

U = r -1 D, (18) 

where T is the covariance matrix with element = CY t (t)Y-(t)- 

To derive the evolution equation for Yj(t;uj), multiply both sides of (10) by Uj(x,i) and take the inner 
product to obtain 

du(x,t) , A / / dYi(t:uj) v—*,,/ \ / duAxA) , ,\ 

y m J ,u 3 (x,t)j^+J2(Mx,t),u ] (x,t)) x ^ J +E r ^")\ dt ,Uj( - X,t) / x (19) 

= (C[u(x,t;uj); 0(w)],Uj(x,t)) x . 

Note that the third term of left hand side in (19) vanishes completely due to the DO condition (11), while, 
the second term vanishes for all i ^ j owing to the orthonormality of Ui(x,t), thus 

^^ + (^^M*,t)) ={C[uix,tiu,)i8{w)],u i (*,t))x- (20) 

5 



Note that multiplying (14) by Ui(x,t) and taking inner product gives 
' du(x, t) 



n/ ,Ui{x,t)J ={E u [C[u(x,t;u]);0(u;)]},u i ( X ,t)} x . (21) 



Using (21) in (20) gives the evolution equation for Yi(t;uj) as 



dt 



(C [u(x, t; to); 0(u)] - E" [C [u(x, t; w); 0(w)]] , Uj(x, t)) x . (22 



3.1.2. Bi- orthogonal Expansion 

Note that the numerical solution of (SPDE) using the DOFE method provide the samples of u(x,t;uj) 
through the coefficients Yi(t;u>), whereas, the Bayesian inference requires analytic form of the probability 
distribution of the system response, thus, the DOFE method can not directly be used for the Bayesian 
inference. In the present paper, a bi-orthogonal expansion approach is proposed to impose the orthogonality 
based geometric structure on the stochastic dimension. Consider a gPC expansion of Yi(i;w) truncated at 
P terms as 

p 

Yi(t;w) =X^(*M»(€M), (23) 
P =i 

where ip p (£(u])) are the orthogonal polynomials from the Askey scheme, while £(w) € L 2 (E) are the random 
variables with appropriate probability density function [6]. Use (23) in (9) to get 

N P 

u(x, t; u) - u(x, t) + Y p Wp(£MHfo *)• (24) 

i— 1 p— 1 

Equation (24) is termed here as the fei-orthogonal expansion. Differentiate (23) with respect to time and 
use the Galerkin projection to obtain 

dYl(t) 1 / \ 
-^T 1 = 77^( (FHx,t;u);0(u)}-E u [F [u(x, t; w); 0(w , u t (x, t)) x , V P «H)) • (25) 
dt W) s / 'a 

Equations (14), (17) and (25) forms dynamically bi-orthogonal field equations (DBFE) that define the 
dynamic evolution of the mean u(x,t), the eigenfield Ui(x,t) and the associated coefficients Yp(i). The 
resultant bi-orthogonal expansion (24) approximates the system response u(x, t\ oj) to an arbitrary accuracy 
depending on the number of eigenfunctions used, N, and the number of expansion coefficients for each 
eigenfunction, P. 

3.1.3. Boundary Conditions 

To define boundary conditions for DBFE, consider a generic Karhunnen-Loeve expansion of h{j3,t\uS) 

N 

h{{3, t; u) = h(J3, t) + Y Y i (*5 w )^(/3' *)• (26) 
i=l 

Applying the expectation operator to (26), boundary condition for the mean is given by 

B(u{x,t)) =h(J3,t). (27) 

By multiplying Yj(t;u>) to (26) and applying the expectation operator to obtain boundary condition for 

Ui(x,t) 

N 

B{Ui(p,t)) ^C-^./ [h^t-^it-u)] . (28) 



3.2. Bayesian Inference 

Without loss of generality, the proposed method is described here for a spatially varying uncertain 
parameter with prior given by a scalar stochastic process v(x;uj), i.e. 0(u) = {v(x;ui)} . Use a KL 
expansion of v(x; to) as 

N 

v(x;uj) =v(x) + ^2^/X l v i (x)xi, (29) 
i=l 

where Xi are independent identically distributed zero- mean random variables, while, A^ and Vi(x) are the 
eigenvalues and eigenfunctions of the covariance function oiv(x;u>). For the covariance function C v (xx,X2), 
Aj and Vi(x) are solution of the eigenvalue problem 

/ C v (x 1 ,x 2 )v i {x 1 )dx 1 = \iVi{x 2 ). (30) 
J x 

For a Gaussian process prior, Xi are standard normal random variables, whereas, for a generic stochastic 
process prior, Xi are given by 

Xi = —7= (v(x; u) - v(x)) Vi{x)dx. (31) 

V *i JX 



Use the gPC expansion of Xi 



Xi = £x^p«M)> (32) 
P =l 



where Xp are the gPC expansion coefficients, in (29) to get the bi-orthogonal expansion of v(x;w) as 

N P 

v(x; uj) = v(x) + V^Vi(n)x],'ipp(£(u))- (33) 

i— 1 p— 1 

The expansion (33) is used in C [u(x, t; lu); 8(lj)] to define the RHS of the DBFE governing equations ((14), 
(17) and (25)). The numerical solution of the resultant DBFE governing equations give the bi-orthogonal 
expansion (24) of the system response u(x, t;u)). 

Remark 2. The numerical solution is initiated with the initial condition for the mean u{x,t) given by 

u(x, t) = E u [F[u(x, 0; w); 0(w)]] , (34) 

while, the initial conditions for the eigenfield are given by 

Ui(x,t) = Vi(x). (35) 

Since the stochasticity in (SPDE) emanates due to the uncertainty in v(x;uS), (SPDE) is initialized with 
a deterministic initial condition. Thus, the initial condition for the expansion coefficients Y* (t) is given by 

r;(0)=0; Vi = l,..,N; p=l,...,P. (36) 
The bi-orthogonal expansion of u(x,t;ui) is used in (8) to define the likelihood 

f(y e | e,of,A,) oc| E \-i cxp Li^ E -i^ ( (37) 

where 

r]= \y e (x k ,t)- [ u(x k ,t) + > " y *Y2(t)ui(x k ,t)ipp(£(u))) I ; k=l,...,M}. (Ss) 



N P \ "> 

l(x k ,t) +^^Y*(t)ui(x k ,tW P (e(w))\ ; k = l,...,M 

i=l p=l / j 



Note that conditional on the hyper-parameters of the discrepancy function, a\ and Ai, £ are the only 
uncertain parameters in (37). Thus, the Bayesian calibration problem is reformulated in the space £ 2 (S) 
as the inference of In the present paper, the proposed method is demonstrated for Hermite polynomials 
as gPC basis, where £ are the independent identically distributed standard normal random variables, thus, 
the prior for £ is given by 

/(^ocffexpf-f), (39) 
fe=i ^ ' 

where N z is the dimension of stochasticity. Using (39) and (37) in (8), the proposed formulation for the 
Bayesian inference is 

malXi) oc| £ |"i exp f-^E^ x Orf)"^" 1 exp (-^) x n^) QAi_1 ex P 



i=l 



Si 



X 

fe=l 



n^p 9 ■ (4o) 



Note that (40) does not involve the solution of T(x,t,0(u>)), thus, the posterior distribution can can be 
explored efficiently using MCMC. In the present paper, Metropolis-Hastings algorithm [14, 31] is used to 
sample from the posterior distribution. 



4. Numerical Example: 2D Transient Diffusion Equation 

The efficacy and efficiency of the proposed method is investigated for calibration of a two-dimensional 
transient diffusion simulator with uncertain source location and diffusivity field. Present paper considers 
a stochastic transient diffusion equation defined over a two-dimensional domain X = [—1,1] x [—1,1] with 
adiabatic boundary conditions as 

g t ' = V (y(x; u)Vu(x, t; u)) + ^ ^ exp f - v L-g > J (41) 

Vu{x,t;u) ■ h = (42) 
u(x,t;u) = 0, (43) 

where v{x]uS) is spatially varying diffusivity field, while total N s source term are active at locations Zi 
with source strength Si. The diffusivity field, v(x;u)), and the location of the source, 2, are assumed to be 
uncertain. Efficacy of the proposed method is demonstrated using the 'hypothetical test bed' data, which is 
defined using the numerical solution of (43) for completely known source location and the diffusivity field. 
In the present paper, the proposed method is demonstrated for a single source located at (0.2, —0.2), which 
is active during time [0, 0.01]. The spatial variation of the diffusivity is assumed to take the form 

v{x) = 0.05(^o + 10.0 + 0.25.x + 0.65y + x 3 + y 3 ), (44) 

where uq is a user defined constant. Figure 1(a) shows spatial variation of the diffusivity. The deterministic 
numerical solution is obtained using a second order central difference scheme in spatial dimension with uni- 
form grid spacing h — 0.02, while the explicit fourth order Runge-Kutta scheme is used for time integration 
with the time step At — 0.0001. Figure 1(b) shows the numerical solution at t = 0.05s, while the solution 
at t = 0.1s is shown in Figure 1(c). Note that the source has peak strength at z; and reduces exponentially 
with the distance, resulting in the peak value of u at the source location and the subsequent diffusion with 
time to other locations. Upon removal of the source, diffusion of u is non-uniform owing to the non-linear 
diffusivity. 
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Diffusivity 



time=0.05s 



time=0.1s 






(a) spatial variation of diffusivity 



(b) u-ficld at t = 0.05s 



(c) u-field at t = 0.1s 



Figure 1: Solution of two-dimensional transient diffusion equation 



4-1- DBFE Formulation 

For notational convenience, define 



S(x;w) 



27TCT 



s I (z — xY 

exp 



2a 2 



(45) 



which is uncertain owing to the uncertainty in the source location z. The prior uncertainty in z is expanded 
in a gPC basis, while the Galerkin projection is used to obtain the resultant gPC coefficients, S(x), of 



P =i 



(46) 



The prior uncertainty in v(x;lo) is represented using a Gaussian process, which is spectrally represented 
using the bi-orthogonal expansion as 



N P 



i—1 p—1 



(47) 



where V(x;uj) is the mean, Vi(x) are the eigenfunctions of the covariance function of v{x;ui) and V£ are the 
respective expansion coefficients. Use (46) and (47) in (43) to obtain the differential operator in (SPDE) 
as 

N P 

C [u(x, t; u); 9(cj)} = V[V(x)Vu(x, t) + V(x) £ ^ ^(t)V» P (€(w))V«i(x, t) 

i—1 p—1 

N N P P 

+ E E E E v} Y i(tM*)Mt(u))Mt(u))v^*, *) (48) 

i—1 j—1 p—1 q—1 
N P P 

+ E E V; Vl (x)^ p {i{uj))Vu(x, t)} + 5(a0lW€M)- 

i= 1 p—1 P—1 

Use (48) in (14), (17) and (25) to obtain the DBFE governing equations for the two-dimensional transient- 
diffusion equation. 
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Figure 2: Comparison of accuracy and computational efficiency of DBFE with gPC and Monte Carlo method 



4-2. Solution of Forward Problem 

The prior uncertainty in the source location is specified using independent Gaussian distribution for x 
and y co-ordinate with mean and standard deviation 0.3, while, the prior for diffusivity v(x\ uS) is specified 
using a Gaussian process with mean 

v{x) = 0.05(^o + 10.0 + 0.25a; + 0.65y) (49) 

and the squared exponential covariance function 

C{x 1: x 2 ) = al exp (-Ai(n - x 2 f - A 2 (j/i - y 2 ) 2 ) , (50) 

where is the variance of the Gaussian process and A, is the correlation length. 

Efficacy and the computational efficiency of the proposed Bayesian inference depends on the ability of 
the DBFE method to accurately solve the forward forward problem. In the present subsection, accuracy and 
the computational cost for the numerical implementation of DBFE is compared against the Monte Carlo 
and the generalized Polynomial Chaos method (see [33] for the gPC formulation of (43)). 

Figure 2 shows the accuracy and computational efficiency of the DBFE and the gPC method for different 
values of number of eigenfunctions used, N, and the order of the polynomial chaos basis, p. The accuracy is 
compared using the Monte Carlo method with 10000 samples, which are collected at the computational cost 
of 6616.17 seconds. The computational cost for solution of the forward problem increases with increase in N 
and p for both the DBFE and gPC method. Note that the stochastic dimension for the present problem is 
N + 2 (N dimensions representing the truncated KL expansion, while 2 dimensions representing uncertainty 
in the source location) , for which the number of polynomial chaos terms is given by 

P (7V + 2)!p! +L (51) 

Since implementation of the gPC method requires numerical solution of P PDEs (see [34] for details), 
whereas, the DBFE method involves numerical solution of (N + 1) PDEs and N x P ODEs, the increase in 
computational cost with p is significantly higher for the gPC method as compared with the DBFE method. 
The computational cost of the gPC method is comparable to the Monte Carlo method for N — 6 and 
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second order polynomial chaos basis, while the computational cost is higher than the Monte Carlo method 
for third order polynomial chaos basis with N > 3, rendering the gPC method computationally intractable. 
The DBFE method is numerically implemented at a computational cost comparable to the gPC method for 
p = 1, while, the computational cost for the DBFE method for p > 2 is considerably lower than the gPC 
method. The ii-error in variance is shown in Figure 2(b). The error is comparable for both the DBFE and 
the gPC method, which decrease with N and reaches the limiting value for N > 4, though the limiting value 
is higher for p = 1. Note that the transient diffusion equation involves multiplication of the diffusivity v 
with Vu, thus, appropriate spectral representation requires use of the second order polynomial chaos basis. 
From the results, it may be concluded that (43) can be numerically solved using the DBFE method at 
significantly lower computational cost than the gPC method with the comparable accuracy. 

4-3. Solution of Source Inversion Problem 

The proposed method is used for inference of the source location and the diffusivity. Prior uncertainty 
in the source location is given by independent Gaussian processes in x and y directions, with Af(0. 0,0.3) 
prior. Prior uncertainty in the diffusivity is specified using the Gaussian process with the mean (49) and 
the covariance function (50) with a 2 = 0.3 and A = 1.5. The prior uncertainty is propagated to the 
system response using the DBFE method with N = 5 and p = 2. The deterministic solution of the 2D 
transient diffusion equation at time t = 0.02 seconds with the source removed at t = 0.01 seconds is used as 
experimental observations. The source is assumed to be located at [0.2,-0.2]. Total 25 uniformly spaced 
data points are used for the Bayesian inference. 1% experimental uncertainty is assumed in each data point. 
The model structure uncertainty is defined by specifying the prior probability distribution for a 2 and X$. 
Inverse Gamma distribution /G(6.0, 2.0) is used for er|, while, the prior for \$ is given by the Gamma 
distribution G(6. 0,2.0). Posterior distribution is explored by collected 10000 samples using MCMC, after 
burnout period of 1000 samples. 

Figure 3 shows the posterior probability density contours for the source location. The contour is shown 
for the posterior density obtained using the proposed method (dashed line) and the direct MCMC sampling 
from the posterior distribution (solid lines). The source location is predicted accurately, while, the posterior 
density obtained using the proposed method agrees closely with the direct MCMC sampling, demonstrating 
the efficacy of the DBFE based Bayesian inference. 

Figure 4(a) shows the Li-error in posterior variance of the diffusivity between the DBFE the direct 
MCMC sampling. The maximum error is of the order of 10~ 3 , indicating the close agreement in variance 
for the posterior probability of the diffusivity obtained using the DBFE and the direct MCMC sampling 
methods. Figure 4(b) shows the Li-error for posterior mean of the diffusivity for the DBFE method obtained 
against the 'true' diffusivity. Note that the error reduces non-uniformly, indicating the effect of the location 
of the experimental observations on the proposed Bayesian inference. 

Figure 5 shows the comparison of the posterior probability distribution of a 2 and Xg obtained using the 
direct MCMC sampling and the proposed method. The posterior distribution for Xg obtained using the 
proposed method matches closely with the direct sampling, however, the match is comparatively poor for 
the posterior distribution of a 2 . Note that the bi-orthogonal expansion obtained using the DBFE method 
acts as an emulator of the 2D transient diffusion equation, which is used in the Bayesian inference as against 
the simulator in the direct MCMC sampling. Thus, any remnant error in the bi-orthogonal expansion is 
considered as the uncertainty in the model structure, resulting in the difference in the posterior probability 
distribution for a 2 . The posterior probability for X$ has moved towards the left for both the cases, indicating 
increased correlation for the model structure uncertainty, while, the posterior distribution for a 2 moves 
towards right, indicating higher posterior confidence on the simulator. 

Figure 6 shows the Li-error in the posterior mean for the system response u, which is defined against the 
true spatial distribution of u. The posterior mean of the system response is calculated by substituting the 
mean of £ in the bi-orthogonal expansion. The maximum Li-error is of the order of 10 _1 , which is located 
in the boundary region where experimental data is not provided for the Bayesian inference. In the region 
where experimental observations are available, the error is significantly low with the minimum value of the 
order of 10~ 7 . 
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Figure 3: Posterior Probability Contours for the Source Location 



5. Concluding Remarks 

The paper has presented a dynamic bi-orthogonality based approach for computationally efficient im- 
plementation of the Bayesian inference. The proposed method can be applied for calibration of a simulator 
represented using partial differential equation with high dimensional uncertainty. Though the method re- 
quires reformulation of the governing equations, existing schemes can be extended in a straightforward 
manner for numerical solution of the DBFE. 

Efficacy of the proposed method is demonstrated for calibration of a two-dimensional transient diffusion 
equation with uncertain source location and the diffusivity. Computational cost of the proposed method 
for uncertainty propagation is compared against the gPC and the Monte Carlo method. Note that for 
low dimensional uncertainty, computational cost of the gPC method is comparable to the DBFE, however, 
as dimensionality of the uncertainty increases, the DBFE method provide the solution of the SPDE at a 
significantly less computational cost than the gPC method with comparable accuracy. Accuracy of the 
proposed method to infer the uncertain parameters is compared against the direct MCMC sampling. The 
method provide accurate inference of the source location with the marginal posterior distribution matching 
closely with the MCMC sampling. The method is found to accurately infer the posterior distribution of the 
spatially/temporally varying parameters. In the present paper, the proposed method is demonstrated in the 
Gaussian context. In the future, efficacy of the proposed method will be demonstrated for a more generic 
non-Gaussian non-stationary setting. 
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Figure 6: Li-error in posterior mean for u 
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